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ABSTRACT 


The  steady  state  problem  of  the  planar  hydrostatic  foil  bearing  is 
analyzed  and  solved  numerically.  Two  techniques  of  solution  are  used.  One 
method  is  simulation  in  time  with  asymptotic  approach  to  steady  state.  This  is 
achieved  by  a preprocessor  which  automatically  sets  up  the  numerical  computer 
program.  The  second  method  is  an  iterative  shooting  technique.  The  results  agree 
well  with  one  another.  Curves  of  pressure  and  typical  film  thickness  versus  flow 
are  presented. 
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NOMENCLATURE 
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G constant  found  from  the  solution  of  Eqn.  (25) 
h local  film  thickness 

F dimensionless  film  thickness  h/r^ 

hg  dimensionless  distance  of  foil  asymptote  to  cylinder 
H stretched  film  thickness  coordinate 

H,P  auxiliary  normalized  values  of  H for  computational  purposes 
m,  n exponents,  to  be  determined 
p local  pressure  in  film 

Pg  ambient  pressure 

Q volume  flow  rate  per  unit  bearing  width 
R local  radius  of  curvature 

r^  cylinder  radius 

t time 

T foil  tension  per  unit  width 

e dimensionless  flow  rate 

•■oT 

H gas  viscosity 

n dimensionless  pressure 

d angular  coordinate  (origin  at  point  of  tangency) 

r dimensionless  time 

( stretched  angular  coordinate 

—A 

auxiliary  normalized  values  of  ( for  computational  purposes 

Subscripts 

c pertains  to  cylinder  center 

g pertains  to  source  location  (groove) 

t pertains  to  point  of  tangency 

oo  pertains  to  end  of  lubrication  zone 
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INTRODUCTION 


Despite  considerable  effort  in  the  field  of  foil  bearings  over  the  past 
two  decades,  an  analysis  of  hydrostatic  foil  bearings  has  not  been  published.  This 
problem  is  of  interest  in  magnetic  recording  and  elsewhere.  The  purpose  of  this 
report  is  to  fill  this  gap. 


The  initial  goal  of  this  investigation  was,  actually,  to  illustrate  the 
applicability  of  a preprocessor  computer  program  for  the  automatic  solution  of 
partial  differential  equations.  Ironically,  it  turned  out,  the  automatic  solution  was 
only  partly  successful  and  numerical  difficulties  were  encountered  in  certain  param- 
eter ranges,  requiring  some  human  intervention.  Consequently,  an  alternative  approach 
was  utilized  and  a complete  set  of  data  curves  generated.  It  is  felt,  however,  that 
the  numerical  experience  gained  in  this  effort  as  well  as  the  concepts  of  the  pre- 
processor  are  of  interest  and,  therefore,  they  are  reported  here  as  well  as  the  solution 
of  the  particular  problem  at  hand. 


HYDROSTATIC  FOIL  BEARING  MODEL 


Let  the  configuration  shown  in  Figure  1 be  studied.  Here,  an  infinitely 
wide,  hydrostatic,  cylindrical,  perfectly  flexible  foil  bearing  with  incompressible 
lubricant  is  depicted  schematically.  In  this  configuration,  one  line  feed  source  is 
located  at  0 = 0g.  A second  line  source  is  symmetrically  situated  at  the  other  side 
of  the  bearing.  It  is  assumed  that  the  feed  pressure  Pg  is  prescribed  and  maintained 
constant.  The  film  thickness  and  pressure  are,  then,  governed  by  the  following 
differential  equations: 


Using  the  dimensionless  representation 


• -«c  (CENTER) 


(SOURCE) 


(TANGENCY) 


STATIONARY 

FOIL 


( END  OF 
LUBRICATION 
REGION  ) 


Figure  1 Cylindrical  Hydrostatic  Foil  Bearing 
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and  restricting  the  problem  to  steady  state,  Eqns  (1),  (2)  may  be  rewritten  as 

_«t£ 

90  V-  -r 


7t 


I-  ... 


(3) 


(4) 


where  Q is  an  integration  constant  representing  the  volume  flow  rate  per  unit  width. 
Eqns.  (3),  (4)  may  be  combined  to: 

% 

& + 


fn 


(5) 


where  the  dimensionless  flow  rate  e is  defined  by 

n,  r- 

It  may  be  observed  that  under  steady  conditions  there  is  no  flow  in  the  region 
< 0 < 0g  (due  to  symmetry  we  restrict  our  treatment  to  one  half  of  the 
problem).  The  pressure  in  this  region  is,  therefore,  uniformly  P=Pg-  In  dimen- 
sionless form: 


TT-TT, 


=jiA 


TM 


The  boundary  conditions  are: 

At  B = e„ 


z -> 

h,  h and  II  {>re  continuous 


(6) 


A, 


Ct«0 


- I 


(7) 

(8) 


Condition  (8)  specifies  that  the  foil  approaches  an  asymptote  located  at  some  small 
unknown  dimensionless  distance  hg  away  from  the  point  of  tangency.  This  unknown 
may  be  eliminated  from  Eqn.  (8)  uy  differentiation. 

Consider  now  a hypothetical  process  in  which  the  tension  is  gradually 
increased  while  0^,  0g,  Pg-Pa*  p remain  unchanged.  Once  anticipates  a corresponding 
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reduction  of  flow,  or  of  e.  (Note  that  this  process  implies  either  a change  in 
restriction  or  in  the  source  pressure,  to  maintain  a constant  Pg  despite  changes 
in  Q).  Following  a derivation  analogous  to  Ref  [1]  one  assumes 

€ 


The  following  requirements  may  be  imposed:  The  two  sides  of  Eqn 
(5)  must  balance  as  e o;  hence. 


4»f  -3>m  s / 


Secondly,  at  least  one  variable  term  in  Eqn  (2)  must  not  vanish  as  e o;  hence, 
7t  - ZHT 


It  follows  then,  that 


'/t 

yt  = Vs- 

Thus,  in  the  region 
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where  is  the  unknown  dimensionless  film  thickness  at  ^ 

The  formulation  of  the  problem  in  region  dg  < 0 becomes 


(17) 


(18) 


(19) 


(20) 


(21) 


(22) 


The  fact  that  we  have  four  boundary  conditions  for  a third  order 
equation  indicates  that  one  of  the  parameters  of  the  problem  is  not  free.  The 
following  functional  form  is,  therefore,  deduced; 


(23) 

(24) 
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RESULTS  AND  DISCUSSION 


The  results  are  summarized  in  Figures  2-5.  They  are  expressed  in 
terms  of  two  dimensionless  parameters,  namely:  ^g/^c  sp«cifving  the  source 

describing  flow  rate.  Figure  2 illustrates  some 

typical  tape  contours  while  Figures  3-5  present  the  dimensionless  groove  pressure 
and  some  characteristic  film  thicknesses. 

The  basic  behavior  of  the  bearing  as  it  appears  from  these  graphs 
may  be  described  by  means  of  a conceptual  experiment.  In  this  experiment, 

Pg-Pg  is  slowly  increased  while  the  geometry  and  tension  remain  constant.  Initially, 
when  the  pressure  is  low,  there  is  no  flow  (Q=0).  When  Pg-Pg  is  increased  just 
above  the  threshold  level  of  T/r^,  the  tape  lifts  off  and  starts  forming  into  a 
"bubble"  shape  (Figure  2).  This  contour  is  typified  by  a large  ratio  of  h^/h^, 
but  the  vanishingly  small  values  of  both  h^  and  h^  make  the  radius  of  curvature 
over  the  source  equal  to  r^.  Increase  in  pressure  above  the  threshold  results  in 
growth  of  h^  as  well  as  h^  and  formation  of  smaller  radius  of  curvature  over  the 
source.  The  rate  of  growth  of  h^  is  slower  than  that  of  h^,  resulting  in  a gradual 
loss  of  the  "bubble"  shape.  Hence,  the  radius  of  curvature  over  the  groove  reaches 
a minimum  and  then  starts  growing, eventually  exceeding  r^.  This  growth  causes 
the  impedance  of  the  bearing  to  decrease  and  the  flow  to  increase.  The  greater 
radius  of  curvature  implies,  however,  that,  at  the  same  time,  the  pressure  required 
to  support  the  foil  tension  is  smaller.  The  above  behavior  results  in  the  unexpected 
phenomenon  that  the  same  value  of  Pg  generates  two  distinct  flow  rates.  This 
should  not  be  surprising,  however,  if  one  remembers  that  the  shape  of  the  flow 
channel,  fornted  by  the  foil,  is  markedly  different  in  the  two  cases. 
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A natural  question  which  arises,  is  what  flow  rate  will  the  physical 
system  select  in  view  of  the  ambiguity.  The  analysis  in  this  work,  determines 
only  the  static  equilibrium  conditions  and  not  their  stability  and,  therefore,  cannot 
provide  a complete  answer  to  the  question.  Nevertheiess,  a limited  insight  will  be 
provided  in  the  following  discussion. 

At  the  onset  let  it  be  noted  that  the  assumption  has  been  implicity 
made,  that  Pg  may  be  controlled  independently  of  Q.  In  any  practical  situation 
the  restrictor  characteristic  presents  a curve  of  pg  vs.  Q,  which  may  be  superimposed 
on  Figure  3.  The  intersection  determines  possible  operating  points  which  are 
schematically  illustrated  in  Figure  6. 

Neglecting  the  effect  of  damping  which  may  be  brought  out  only  by  a 
dynamic  analysis,  it  may  be  shown  that  intersections  of  type  A,  C tend  to  be  stable 
whereas  intersections  of  type  B tend  to  be  unstable.  The  argument  is  essentially  as 
follows.  Considering  point  B for  example,  a downward  fluctuation  in  Pg  causes 
the  bearing  outflow  to  exceed  the  supply  through  the  restrictor.  This  tends  to 
push  the  operating  point  further  down  from  B by  causing  a deficiency  in  fluid 
at  the  bearing  inlet.  Thus,  for  an  intersection  of  restrictor  characteristic  and  film 
characteristic  such  as  that  shown  at  B,  a disturbance  in  Pg  away  from  the  equili- 
brium point  results  in  a further  excursion  of  the  operating  point  in  the  same 
direction.  However,  for  intersections  of  types  A or  C,  disturbances  in  Pg  from 
the  equilibrium  point  result  in  a restorative  effect,  suggesting  stability. 

Another  observation  of  interest  is  that  » 0.05  is  a threshold  value, 
below  which  contact  between  the  tape  and  the  cylinder  may  arise.  In  these  situations 
the  source  is  located  too  close  to  the  point  of  tangency  for  useful  operation.  When 
the  flow  is  minute,  a noncontact  situation  is  possible.  For  slightly  higher  flow  rates, 
the  seal  effect  at  the  tangency  point  disappears  and  flotation  of  the  foil  cannot 
be  maintained  throughout  the  wrap  angle. 
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APPENDIX  I:  PREPROCESSOR  SOLUTION 


PURPOSE  OF  PREPROCESSOR 


With  present-day  powerful  computers,  the  engineer  has  at  his  disposal 
an  invaluable  tool  for  optimization  and  evaluation  of  design  parameters.  Frequently, 
this  evaluation  is  based  on  numerical  simulation  of  the  behavior  of  a physical  model 
represented  by  a partial  differential  equation  (POE).  In  order  to  perform  the  simulation, 
it  becomes  necessary,  in  such  cases,  to  discretize  the  partial  differential  equation  (POE) 
and  to  write  a computer  program  that  will  solve  the  resulting  difference  equation 
and  print  or  display  the  results  numerically  and/or  graphically.  This  task  is,  often, 
time-consuming  and  usually  requires  several  debugging  runs. 

Though  laborious,  this  programming  effort,  commonly,  consists  of 
rather  similar  subtasks.  In  general,  the  programming  job  may  be  described  as  the 
insertion  of  problem  dependent  details  into  a basic  program  template  designed  for 
the  particular  algorithm  and  a broad  class  of  equations.  Consequently,  it  is  feasible 
to  carry  out  this  process,  too,  with  the  aid  of  the  computer.  This  is  done  by 
means  of  a preprocessor  or  a precompiler.  The  precompiler  is  a computer  program 
which  accepts  an  input  describing  the  partial  differential  equation  and  its  boundary 
conditions  in  a concise  and  '-eadable  mathematical  symbolism,  and  outputs  a high 
level  language  (e.g.,  PLI)  program  which  is  then  compiled  by  the  language  compiler 
just  like  a human-written  program.  The  compiler-generated  object  program  is 
finally  executed.  (Figure  A1). 


■Nai  Ai'  Wlr,it&. 
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( Program  y 
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Results 


Figure  A1  Schematic  View  of  the  Role  of  the  Precompiler  in  Problem  Solution. 
(Broken  lines  indicate  conventional  computer  solution.  Solid  lines 
add  the  effect  of  the  precompiler.) 


In  this  report,  we  would  limit  the  discussion  to  a precompiler  for  the 
class  of  parabolic,  high  order,  one  dimensional  partial  differential  equations. 

We  will  restrict  ourselves,  further,  to  the  implicit  solution  algorithm.  This  class  of 
problems  is  still  rather  broad;  and  in  particular,  it  is  applicable  to  many  foil- 
bearing problems. 
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FORMULATION 


In  order  to  allow  solution  by  means  of  a program  generated  by  our 
precompiler  for  parabolic  partial  differential  equations,  the  formulation  (18)  - (22) 
for  H and  Ilg  will  be  replaced  by  the  following  alternative  formulation  defining  a 
function  H (|,r).  H will  later  be  transformed  to  the  desired  function  H ({). 


-a>V  ^ 


f.r# 


When  a steady  state  solution,  ^ 3 ^ J 

satisfying  these  requirements  is  obtained,  a parameter  G,  elaborated  upon  later,  may 


be  evaluated  from 


(30) 


AMPBX 


^ and  may  then  be  obtained  from  the  numerical  solution 
obtained  by  the  transformations 


//.  w/<r 


(31) 


(32) 


3V 


7 • ,T 


(33) 


Substitution  of  the  transformations  (31)  - (34)  into  Eqns  (25)  - (29)  will  verify 
that  the  function  H indeed  satisfies  all  the  requirements  (18)  - (22).  The  equality 
of  the  ratios  ^ ^ simplifies  the  presentation  of  the 

results. 


(34) 


The  value  of  G,  though  theoretically  a constant,  may  be  expected  to 
vary  somewhat  with  ( due  to  truncation  and  round-off  errors  in  the  solution.  It 
was  rather  unexpected,  however,  to  find  that  in  certain  parameter  ranges  unacceptable 
nonuniformities  in  G were  found  numerically.  However,  for  those  parameter  values 
for  which  the  values  of  G were  sufficiently  uniform,  the  results  agreed  very  well  with 
those  of  the  alternative  technique  described  below. 


AMPEX 


PRECOMPILER  INPUT 


The  precompiler  input  corresponding  to  Eqns.  (25)  - (29)  is  shown 
in  Figure  A2.  The  first  line  identifies  the  user  supplied  program  name  EXTP, 

This  name  is  used  as  a prefix  in  naming  subroutines  generated  by  the  preprocessor 
to  facilitate  readability  of  the  generated  program  by  the  user.  Next,  room  is 
provided  for  some  options.  In  the  present  precompiler  version,  this  is  superfluous 
since  only  one  option  now  exists.  Next,  one  may  note  in  Figure  A2  some  text 
enclosed  between  /*  . . . */.  This  text  constitutes  a comment  which  is 

ignored  by  the  precompiler.  Following  this  text,  some  declarations  appear.  These 
declarations  categorize  identifiers  as  parameters,  independent,  dependent  and  entry 
names.  Parameters  signify  variables  whose  numerical  value  will  be  read  in  at  run 
time  by  a standard  read  statement.  Entry  names,  not  illustrated  in  Figure  A2,  are 
names  of  library  or  user  defined  subroutines,  and  so  on. 

Following  the  preamble,  the  actual  statements  of  the  differential 
equation  (DE),  boundary  conditions  (BC),  and  initial  conditions  (1C)  are  provided^ 

The  specification:  BC  LOW  (XI R * XIC,  *)  means:  This  is  a boundary  condition 
at  the  low  end  of  the  interval  of  independent  variable  #1;  i.e.,  X whose  value  at 
that  point  is  the  expression  XI R * XlCt  In  addition,  the  BC  is  valid  for  any 
value  (*)  of  independent  variable  #2,  i.e.,  T.  The  rest  of  the  input  is  self-explanatory. 


t Naturally,  this  value  could  have  been  a constant  or  a variable  name. 
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EXTP:  PROCEDURE  OPTI ONS( PARABOLl C) ; 

I*  EXTERNALLY  PRESSURIZED  FOIL  BEARING  ANALYSIS. 
STATIONARY,  INCOMPRESSIBLE,  PLANAR,  PERFECTLY 
FLEXIBLE,  REFERENCED  TO  A CYLINDRICAL  SURFACE. 

TWO  SOURCES,  SYMMETRICALLY  LOCATED  ABOUT  CENTER.  */ 


DECLARE 

XIC 

PARAMETER, 

XIR 

PARAMETER, 

XI8 

PARAMETER, 

HG 

PARAMETER, 

HMIN 

PARAMETER; 

X 

T 

H 


/* 

/* 
/* 
/* 
/* 

DECLARE 

INDEPENDENT(l), 

INDEPENDENT(2), 


COORDINATE  OF  CENTER  */ 

XIG/XIC  , WHERE  X I G=COORD I NATE 
COORDINATE,  END  OF  LUBRICATION 
H VALUE  AT  SOURCE  */ 

INITIAL  VALUE  OF  M AT  TANGENCY  POINT  */ 


OF  SOURCE  *! 
ZONE  */ 


/*  1ST  VARIABLE  IN  LIST,  SPATIAL  COORDINATE 
/*  2ND  VARIABLE  IN  LIST,  TIME  COORDINATE  */ 


*/ 


DEPENDENT;  /*  DIMENSIONLESS  FILM  THIKMESS,  H_BAR  IN  PAPER  */ 


DE: 

BC  LOW  (XIR*XIC,*); 
BC  LOW  (XIR*XIC,*): 
BC  HIGH  (XI8,*); 

BC  illGH  (XI8,*); 

1C: 

^ END; 


04X(H)*(-H**3)  +@3X(H)*(-3.*H**2*@X(H))*@T(H); 
H«HG; 

02X(H)*(XIR-1. )*XIC+0X(H)*(-1. )=0.; 

0X(H)=XI8; 

02X(H)*1.; 

H»HMlN+X**2/2.; 


Figure  A2  Precompiler  Input  (Corresponding  to  Eqns  (25)-|29)) 
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PRECOMPILER  STRUCTURE 


The  precompiler  is  composed  of  three  sub-programs:  1.  Lexical, 

analyzer,  2.  Syntactic  and  semantic  analyzer;  and  3.  Synthesizer. 

The  lexical  analyzer  consists  of  a small  driver  section  which  picks 
an  input  character,  decodes  it,  and  branches  to  several  possible  routines  accordingly. 
These  routines  identify  literals  (e.g.,  2.0,  20.E-1,  2),  identifiers  (e.g.,  XI R, 

SORT),  reserved  words  (e.g.,  DE,  BC),  comments,  etc.  Identifiers  and  literals 
are  stored  in  a table.  The  output  of  the  lexical  analyzer  is  a uniform  stream  of 
pairs  of  descriptors,  the  first  of  which  indicates  the  type  of  element  encountered 
(e.g.,  **,  +,  identifier).  The  second  descriptor  points,  when  necessary,  to  a table 
of  additional  information  (e.g.,  the  actual  identifier's  name). 

The  syntactic/semantic  analyzer  scans  the  pairs  emitted  by  the 
lexical  analyzer  and  parses  them.  Parsing  means,  in  principle,  matching  of  a 
precompiler  input  statement  to  a set  of  prototype  statements  describing  the 
grammar  of  the  input  language.  For  example,  the  two  prototype  statements 

assignment  <=  variable  = expression 

expression  term  [+  term] 

where  the  brackets  indicate  zero  or  more  repetitions  of  elements  of  the  enclosed 
form,  may  be  matched  to  a fortran  statement  such  as 

A = B + C 

Several  parsing  methods  are  available  in  compiler  design.  Parsing  was  done  in  this 
work  by  recursive  descent  [2j.  In  essence,  for  each  of  the  grammar  nonterminal 
symbols  (e.g.,  assignment,  variable,  term), there  is  a logical  recursive  subroutine  which 
returns  "true"  or  "false"  depending  on  whether  the  terminal  language  symbol 
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(e.g..  A,  B,=, +)  can  be  matched  to  the  prototype.  A subroutine  called  "assign- 
ment" calls  the  subroutine  "variable";  and  if  it  is  true  that  the  symbol  is  a 
variable,  checks  ^whether  the  next  symbol  in  the  input  string  is  a "=",  and  so  on. 
This  parsing  process  is  continued  until  the  input  structure  is  matched  to  the 
syntax  rules  or  else  there  is  no  match  due  to  a user  grammar  error.  Once  a match 
is  achieved,  reference  is  made  to  a set  of  routines  which  produce  the  desired 
output  segment  for  the  final  high-level  language  program.  For  example,  since  the 
differential  equation  is  converted  into  a set  of  algebraic  equations  in  matrix  form, 
the  routines  will  generate  assignment  statements  for  the  coefficients. 

The  synthesis  phase  collects  the  segments  generated  by  the  syntactic/ 
semantic  phase  and  inserts  them  into  a prearranged  program  template.  The  template 
consists  of  a high-level  language  (PLI)  program  for  solution  of  the  desired  class 
of  equations  by  an  implicit  algorithm.  The  template  consists,  basically,  of  the 
following  slots: 

Main  program  ~ a driver  routine  containing  a slot  for  user  input. 

Matrix  generator  - sets  up  the  coefficient  matrix  for  PDE  solution 

and  for  BC. 

Initializer  - the  particlar  initial  conditions  are  set  up. 

Output  - fixed  format  output  of  PDE  solution  is  set  up  with  some 

user  control  over  the  amount  of  printout). 

The  generated  program  uses  a diagonal  band  matrix  solver  which  is  a fixed  sub- 
routine. This  routine  takes  as  its  input  the  matrix  equation  and  its  bandwidth  and 
provides  the  solution. 

With  standard  catalogued  procedures  (such  as  we  have  developed  for 
OS  360),  the  operation  of  the  program  is  quite  simple.  The  input  has  to  be  stored 
in  a data  set.  Data  sets  may  optionally  be  provided  by  the  user  for  storing  the  gen- 
erated PLI  program*,  the  object  program,  and  the  output.  The  precompiler  is 
written  in  PLI  and  is,  thus,  machine  independent  for  all  machines  for  which  PLI 
compiler  is  available. 

*A  listing  of  the  generated  computer  program  is  given  in  Appendix  III. 
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APPENDIX  II:  O.D.E.  SOLUTION 


Let  Eqns  (18)  • (22)  be  reformulated  in  terms  of  alternative  variables 

A A 

H ((),  in  order  to  pernrtit  solution  as  an  ordinary  DE.  Later  the  solution 
for  H will  be  transformed  back  into  H (().  The  reformulated  problem  is: 


u 

0^—  = 1 
” if 

(35) 

A 

At  ^ ® 0 

u 

(36) 

o» 

<X 

II 

(37) 

! 

t 1 
\ 

{ 

A»  A» 
H=Hg 

(38) 

A 

Integration  as  an  initial  value  problem  is  terminated  when  H"  tends  to  approach  a 
constant  say  f)"oo-  The  exact  termination  value  is  not  critical,  but  it  should  be 
sufficiently  large  so  that  Aoo/Hg  » 1.  It  may  be  verified  that  the  conversion 
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transforms  the  solution  aposteriori  to  that  of  the  original  coordinates.  The  parameters 
of  the  problem  which  has  been  solved  are  found  by: 


% M 

A / A 


^ - I H*'"' 


i*  l/C- 


i 
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EXTP: 


APPENDIX  III:  LISTING  OF  COMPUTER  PROGRAM  GENERATED 

BY  PRECOMPILER 


PROCEDURE  OPTIONS  (MAIN); 

DECLARE 

01$  BINARY  FIXED(15,0)INITIAL  (1),/*PRINT  INCR*/ 


0X1$  DEC  FL0AT(6)  I Ml Tl AL( 0. 1 ),/*SPATI AL  STEP*/ 

DT1$  DEC  FL0AT(6)  INITIAL(0. 1),/*TIME  STEP  */ 

IK$  BINARY  FIXED(15,n)  IfMTIAL(l)  ,/*PRINT  FREOUENCY  */ 

IL$  BIN  FIXED(15,0)  I N I Tl AL( 1 ),/*SPATI AL  COUNTR  LO*/ 

IH$  BIN  FIXED(15,0)INITIAL(11),/*SPATIAL  COUNTER  HI*/ 

JFIM  BINARY  FIXED  (15,0)  INITIAL  (1),/*N0  OF  STEPS*/ 

JS$  (KLIN$)  BINARY  FIXED  (15,0),  /*LINE  NUMBER  */ 

KBC$  BINARY  FIXED  (15,0)  INITIAL  ( 2),/*HALF  BAND  OF  MATRIX 

KBCP1$  BINARY  FIXED  (15,0)  INITIAL  ( 3), 

KLIN$  BINARY  FIXED  (15,0)  INITIAL  (0),  /*LINES  STORED*/ 
ORDERS  BINARY  FIXED  (15,0)  INITIAL  ( 4); 

LI:  GET  DATA  ( 

XIC,HG,XIR,HMIN,XI8, 

JFIN, IL$, IH$,DI$,DT1$); 

CALL  HD_EXTP; 

GO  TO  LI; 

HD_EXTP:  PROCEDURE; 

DECLARE 

A$  ( I L$-0RDER$/2: IH$+ORDER$/2,-ORDER$/2:ORDER$/2+l) 

DEC  FL0AT(16),/*MATRIX  EQN*/ 

EMATRXP  ENTRY,  /*SOLVE  MATRIX  EQN  */ 

ESTOREP  ENTRY, /*STORE  FOR  PRINTING*/ 

EWRITEP  ENTRY  ,/*  PRINT  ROUTINE  */ 

H (1:KLIN$, IL$: IH$)  DEC  FL0AT( 16 ) ; /*DEPENDENT  VAR*/ 

A$»0.; 

XL$-XIR*XIC; 

XH$=XI8; 

DX1$-(XH$-XL$)/(IH$-IL$); 

PUT  DATA  ( 

XIC,HG,XIR,HMIN,XI8, 

JFIN, IL$, IH$,DI$,DT1$); 

PUT  DATA  (DX1$); 

I LP$-I L$+KBC$; 

IHM$-IH$-KBC$; 

DX2$-DX1$*DX1$;DX3$-DX2$*DX1$;DX4$-DX3$*DX1$; 

CALL  IC_EXTP; 
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L2:  DO  JR$-1  TO  JFIN  ; 

CALL  EQ_EXTP; 

CALL  BC_EXTP; 

CALL  EMATRXP  ; 

IF  (MOD(JR$, IK$)-n)  THEN  CALL  ESTOREP; 
IF((K$-KLIN$)  I (JR$»JFIN))  THEN  CALL  EWRITEP; 
END  L2; 

RETURN; 

IC_EXTP:  PROCEDURE; 

JR$-0; 

K$-l; 

JS$(K$)-JR$; 

DO  l$-IL$  TO  IH$; 

HCK$J$)-HMIN+(l$-IL$)*DXl$**2/2.; 

END; 

RETURN; 

END  IC_EXTP; 

EQ_EXTP;  PROCEDURE; 

DO  l$-ILP$  TO  IHM$  ; 

C4$-(-H(K$,l$)»*3); 

C3$-(-3.*H(K$^ I$)**2*(+1./DX1$*H(K$; l$+0) 

)); 

A$(l$,-2)=(+l./DX4$*C4$-.5/DX3$*C3$)*OTl$; 

A$(l$,-l)-(-4./DX4$*C4$+l./DX3$*C3$)*DTl$; 

A$( I $^+0)«(+F./DX4$*C4$+n.*C3$)*DTl$-l.; 

A$( l$,+l)»(-4./DX4$*C4$-l./DX3$*C3$)*DTl$; 
A$(l$,+2)=(+l./DX4$*C4$+.5/DX3$*C3$)*DTl$; 
A$(IS,+3)»»H(K$,  1$); 

END; 

RETURN; 

END  EfUEXTP; 

BC_EXTP:  PROCEDURE; 

A$(IL$,+0)-+l.; 

A$(IL$,+l)-+0.; 

A$( I L$,+2)-0.; 

A$(l L$,KBCP1$)-HG; 

C2$-(XIR-1. )*XIC; 

Cl$-(-l. ); 

A$(IL$+l,-l)-+2./DX2$*C2$-1.5/DXl$*Cl$; 

A$(l L$+l,+0)-“ 5./DX2$*C2$+2./DXl$*Cl$; 
A$(IL$+l,+l)»+4./DX2$*C2$-0.5/DXl$*Cl$; 

A$(l  L$  + l,+2)  — l./DX2$*C2$; 

A$(IL$+l^KBCPl$)-0.; 

A$(IH$,-2)-+0.5/DXl$; 

A$(IH$,-D  — 2./DXl$; 

A$(IH$,+0)-+1.5/DXl$; 

A$(IH$^KBCP1$)-XI8; 

A$(IH$-l,-2)— l./DX2$; 

A$(IH$-l,-l)-+4./DX2$; 

A$(IH$-l,+0)  — 5./DX2$; 

A$(IH$-l,+l)-+2./DX2$; 

A$(IH$-1,KBCP1$)-1.; 

RETURN; 

END  BC_EXTP; 
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EMATRXP:  PROCEDURE; 

DECLARE 

(l$,KH$,KV$)  BINARY  FIXED  (15,0);  /*L00P  INDICES  */ 
PSl:  DO  l$»IL$  TO  IH$; 

DO  KH$-1  TO  KBCP1$; 

A$(I$,KH$)«  A$(l$,KHS)/A$(l$,0); 

DO  KV$-1  TO  KBC$; 

IF  (KH$«KBCP1S) 

THEN  A$( I $+KV$, KBCPIS )-A$ ( I $+KV$, KBCP1$ ) 
-A$(I$,KBCP1$)*A$(I$+KV$,-KV$); 

ELSE  A$ ( I $ + KV$, KH$-KV$ )-A$( I $+KV$, KH$-KV$  ) 
-A$ ( I $, KH$ )*A$ ( I $ + KV$, -KV$  ) ; 

END; 

END; 

END  PSl; 

PS2;  DO  l$»IH$  TO  IL$+1  BY  -1; 

DO  KV$«1  TO  KBC$; 

A$( I $-KV$, KBCP1$ )-A$ ( I $-KV$, KBCP1$ ) 
-A$(I$,KBCP1$)*A$(I$-KV$,KV$); 

END; 

END  PS2; 

RETURN; 

END  EMATRXP; 

ESTOREP:  PROCEDURE; 

DECLARE 

1C  BINARY  FIXED  (15,0)  ; 

K$»l+  MOD(K$, KLIMS); 

JS$(K$)-JR$; 

DO  l$-IL$  TO  IH$; 

H(K$, 1$)*  A$(I$,KBCP1$); 

END; 

RETURN; 

END  ESTOREP; 

EWRITEP:  PROCEDURE; 

DECLARE 

(IC,KC)  BINARY  FIXED  (15,0); 

PUT  EDIT  (• l$/JS$*,(JSS(IC)  DO  IC-1  TO  KS ) ) 

(COL(l),A(6),(KS)  F(12,0),  SKIP(2)); 

DO  IC-ILS  TO  IMS  BY  Dl$; 

PUT  EDIT  (IC,(H(KC, 1C)  DO  KC-1  TO  K$ ) ) 

(COL  (1),F(6,0),  (KS)  E(12,4)); 

END; 

PUT  PAGE; 

RETURN; 

END  EWRITEP; 

END  HD_EXTP; 

END  EXTP; 
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